function [cfg, dats]=osl_stats_multiplotER(S)

% [cfg, dats]=osl_stats_multiplotER(S)
%
% calls ft_multiplotER (passing on cfg settings) to display the tstats for contrast S.contrast for
% OAT S.stats
%
% Input:
% S.oat
% S.stats_fname 
% S.modality, e.g. 'MEGPLANAR'
% S.first_level_contrast. First level contrast indices - there can be multiple numbers of these
% S.group_level_contrast. Specify single group contrast index (if group level stats)
% S.cfg 
% S.view_cope. Binary flag, if 0 (default) then tstats are viewed
%
% MWW 2012

global OSLDIR;

try, oat=S.oat; catch, error('No S.oat supplied'); end;
try, stats_fname=S.stats_fname; catch, error('No S.stats_fname supplied'); end;
try, modality=S.modality; catch, modality='MEGPLANAR'; end;
try, contrast=S.first_level_contrast; catch, error('No S.first_level_contrast supplied'); end;
try, group_contrast=S.group_level_contrast; catch, group_contrast=1; end; % specify group contrast if group level stats
try, view_cope=S.view_cope; catch, view_cope=0; end; 
try, cfg=S.cfg; catch, cfg=[]; end;

try, D_tstat=S.D; catch
    S4=[];
    S4.oat=oat;
    S4.stats_fname=stats_fname;
    S4.write_copes=1;
    
    [D_tstat, D_cope]=osl_save_spm_stats(S4); 
    
    stats=osl_load_oat_results(oat,stats_fname);

    if(stats.level==2),
        D_tstat=D_tstat{group_contrast};
        D_cope=D_cope{group_contrast};
    end;
end;

cope=stats.cope;
if(stats.level==2),
    cope = cope(:,:,:,:,group_contrast);
end;

for cc=1:length(contrast),

    data = [];
    data.time = stats.times;
    data.dimord='chan_time';

    if(strcmp(modality,'MEGMAG')),
        cfg.layout      = [OSLDIR '/layouts/neuromag306mag.lay'];  
        chanind = strmatch('MEGMAG', D_tstat.chantype); 

        data.label=D_tstat.chanlabels(chanind);

        ts=zeros(102,length(data.time));
        for ind=1:length(chanind),
            indp=(ind-1)*3+3;
            if(view_cope)
                ts(ind,:)=cope(indp,:,contrast(cc));
            else
                ts(ind,:)=cope(indp,:,contrast(cc))./stats.stdcope(indp,:,contrast(cc));
            end;
        end;

    elseif (strcmp(modality,'MEGPLANAR')),

        cfg.layout      = [OSLDIR '/layouts/neuromag306mag.lay'];  
        chanind = strmatch('MEGMAG', D_tstat.chantype); 

        data.label=D_tstat.chanlabels(chanind);

        ts=zeros(102,length(data.time));
        for ind=1:length(chanind),
            indp=(ind-1)*3+1;
            if(view_cope)
                ts(ind,:)=permute((cope(indp,:,contrast(cc))+cope(indp+1,:,contrast(cc))),[2 1 3]);
            else            
                ts(ind,:)=permute((cope(indp,:,contrast(cc))+cope(indp+1,:,contrast(cc)))./sqrt((stats.stdcope(indp,:,contrast(cc)).^2+stats.stdcope(indp+1,:,contrast(cc)).^2)),[2 1 3]);
            end;
        end;
        warning('The planar grads are being overlayed on the magnetometer locations and using the magnetometer labels');

    elseif (strcmp(modality,'MEG')),

        cfg.layout      = [OSLDIR '/layouts/CTF275.lay'];      
        chanind = strmatch('MEG', D_tstat.chantype); 

        data.label=D_tstat.chanlabels(chanind);

        ts=zeros(length(chanind),length(data.time));
        for ind=1:length(chanind),
            indp=ind;
            if(view_cope)
                ts(ind,:)=cope(indp,:,contrast(cc));
            else
                ts(ind,:)=cope(indp,:,contrast(cc))./stats.stdcope(indp,:,contrast(cc));
            end;
        end;
        
    elseif(strcmp(modality,'EEG')),   % added by DM
        
        B=sensors(D_tstat,'EEG');
        mat=spm2fieldtrip(D_tstat);
        cfgx=[];
        cfgx.elec.pnt=B.elecpos;
        cfgx.elec.label=B.label;
        lay = ft_prepare_layout(cfgx,mat);
        
        cfg.layout =lay;
        chanind = strmatch('EEG', D_tstat.chantype);
        
        data.label=D_tstat.chanlabels(chanind);
        
        ts=zeros(length(chanind),length(data.time));
        for ind=1:length(chanind),
            if(view_cope)
                ts(ind,:)=cope(ind,:,contrast(cc));
            else
                ts(ind,:)=cope(ind,:,contrast(cc))./stats.stdcope(ind,:,contrast(cc));
            end;
        end;
        
    else
        error('Unsupported modality');
    end;
   
    data.avg = ts;
    data.datatype='timelock';

    dats{cc}=data;
end;

if(~isfield(cfg,'interactive')),
    cfg.interactive = 'yes';
end;
if(~isfield(cfg,'xlim')),
    cfg.xlim        = [stats.times(1) stats.times(end)];
end;

cfg.comment     = '';

if(length(dats{1}.time)>1)
    
    figure; ft_multiplotER(cfg,dats{:});
else

    if(length(dats)==1)
        dats=dats{1};
    end;
    figure; ft_topoplotER(cfg,dats);
    
end;

